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Abstract 

Ultrasound modulation of electrical or optical properties of materials offers the 
possibility to devise hybrid imaging techniques that combine the high electrical 
or optical contrast observed in many settings of interest with the high resolution 
of ultrasound. Mathematically, these modalities require that we reconstruct a 
diffusion coefficient cr(x) for x £ X, a bounded domain in M n , from knowledge of 
a{x)\Vu\ 2 {x) for x £ X, where u is the solution to the elliptic equation — V-aVu = 
in X with u = f on dX. 

This inverse problem may be recast as a nonlinear equation, which formally 
takes the form of a O-Laplacian. Whereas p— Laplacians with p > 1 are well- 
studied variational elliptic non-linear equations, p = 1 is a limiting case with a 
convex but not strictly convex functional, and the case p < 1 admits a variational 
formulation with a functional that is not convex. In this paper, we augment the 
equation for the O-Laplacian with full Cauchy data at the domain's boundary, 
which results in a, formally overdeter mined, nonlinear hyperbolic equation. 

The paper presents existence, uniqueness, and stability results for the Cauchy 
problem of the O-Laplacian. In general, the diffusion coefficient cr(x) can be stably 
reconstructed only on a subset of X described as the domain of influence of the 
space- like part of the boundary dX for an appropriate Lorentzian metric. Global 
reconstructions for specific geometries or based on the construction of appropriate 
complex geometric optics solutions are also analyzed. 

1 Introduction 

Electrical Impedance Tomography (EIT) and Optical Tomography (OT) are medical 
imaging modalities that take advantage of the high electrical and optical contrast ex- 
hibited by different tissues, and in particular the high contrast often observed between 
healthy and non-healthy tissues. Electrical potentials and photon densities are mod- 
eled in such applications by a diffusion equation, which is known not to propagate 
singularities, and as a consequence the reconstruction of the diffusion coefficient in such 
modalities often comes with a poor resolution (3J El l3T) . 
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Ultrasound modulations have been proposed as a means to combine the high contrast 
of EIT and OT with the high resolution of ultrasonic waves propagating in an essen- 
tially homogeneous medium [32]. In the setting of EIT, ultrasound modulated electrical 
impedance tomography (UMEIT), also called acousto-electric tomography, has been 
proposed and analyzed in [21 [TJ Q31 [TBI E21 GS]- In the setting of optical tomography, a 
similar model of ultrasound modulated tomography (UMOT), also called acousto-optic 
tomography, has been derived in [10] in the so-called incoherent regime of wave propa- 
gation, while a large physical literature deals with the coherent regime [U [221(32], whose 
mathematical structure is quite different. 

Elliptic forward problem. In the setting considered in this paper, both UMEIT and 
UMOT are inverse problems aiming to reconstruct an unknown coefficient cr(x) from 
knowledge of a functional of the form H(x) = a(x)\'Vu\ 2 (x), where u(x) is the solution 
to the elliptic equation: 



Here, X is an open bounded domain in IR n with spatial dimension n > 2. We denote by 
dX the boundary of X and by f(x) the Dirichlet boundary conditions prescribed in the 
physical experiments. Neumann or more general Robin boundary conditions could be 
analyzed similarly. We assume that the unknown diffusion coefficient a is a real-valued, 
scalar, function defined on X. It is bounded above and below by positive constants 
and assumed to be (sufficiently) smooth. The coefficient u{x) models the electrical 
conductivity in the setting of electrical impedance tomography (EIT) and the diffusion 
coefficient of particles (photons) in the setting of optical tomography (OT). Both EIT 
and OT are high contrast modalities. We focus on the EIT setting here for concreteness 
and refer to a as the conductivity. 

The derivation of such functionals H(x) from physical experiments, following similar 
derivations in [TJ QUI [23], is recalled in section [2] For a derivation based on the focusing 
of acoustic pulses (in the time domain), we refer the reader to [2]. This problem has 
been considered numerically in [21 [161 123] . In those papers, it is shown numerically that 
UMEIT allows for high-resolution reconstructions although typically more information 
than one measurement of the form H(x) = cr(x)\Vu\ 2 (x) is required. 

Following the methodology in [2] where the two dimensional setting is analyzed, 
[7] analyzes the reconstruction of a in UMEIT from multiple measurements at least 
equal to the spatial dimension n. The stability estimates obtained in [7] show that the 
reconstructions in UMEIT are indeed very stable with respect to perturbations of the 
available measurements. Such results are confirmed by the theoretical investigations in 
a linearized setting and the numerical simulations proposed in |23j. In this paper, we 
consider the setting where a unique measurement H(x) = cr(x)\Vu\ 2 (x) is available. 

The inverse problem as a p— Laplacian. Following [21 HOI [16], we recast the inverse 
problem in UMEIT as a nonlinear partial differential equation; see ([7]) below. This 
equation is formally an extension to the case p = of the p— Laplacian elliptic equations 



V ■ a(x)Vu = in X, 



u = f on dX. 



(1) 
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posed on a bounded, smooth, open domain X C M n , n > 2, with prescribed Dirichlet 
conditions, say. When 1 < p < oo, the above problem is known to admit a variational 
formulation with convex functional J[Vw] = J x H(x)\'Vu\ p (x)dx, which admits a unique 
minimizer, in an appropriate functional setting, solution of the above associated Euler- 
Lagrange equation [TB] . 

The case p = 1 is a critical case as the above functional remains convex but not 
strictly convex. Solutions are no longer unique in general. This problem has been 
extensively analyzed in the context of EIT perturbed by magnetic fields (CDII and 
MREIT) [2U [27J [2H], where it is shown that slight modifications of the 1— Laplacian 
admit unique solutions in the setting of interest in MREIT. Of interest for this paper 
is the remark that the reconstruction when p = 1 exhibits some locality, in the sense 
that local perturbations of the source and boundary conditions of the 1— Laplacian do 
not influence the solution on the whole domain X. This behavior is characteristic of a 
transition from an elliptic equation when p > 1 to a hyperbolic equation when p < 1. 

The inverse problem as a hyperbolic nonlinear equation. When p < 1, the 
above functional J[Vw] is no longer convex. When p = 0, it should formally be replaced 
by J[Vit] = J x H(x) In \Vu\(x)dx, whose Euler-Lagrange equation is indeed ^ below. 
The resulting 0— Laplacian is not an elliptic problem. As we mentioned above, it should 
be interpreted as a hyperbolic equation as the derivation of ^ below indicates. 

Information then propagates in a local fashion. Moreover, compatible boundary 
conditions need to be imposed in order for the hyperbolic equation to be well-posed 
[2"T| |2"9] . We thus augment the nonlinear equation with Cauchy boundary measurements, 
i.e., Dirichlet and Neumann boundary conditions simultaneously. As we shall see in the 
derivation of UMEIT in the next section, imposing such boundary conditions essentially 
amounts to assuming that a(x) is known at the domain's boundary. This results in an 
overdetermined problem in the same sense as a wave equation with Cauchy data at 
time t = and at time t = T > is overdetermined. Existence results are therefore 
only available in a local sense. We are primarily interested in showing a uniqueness 
(injectivity) result, which states that at most one coefficient a is compatible with a 
given set of measurements, and a stability result, which characterizes how errors in 
measurements translate into errors in reconstructions. Redundant measurements then 
clearly help in such analyses. 

Space-like versus time-like boundary subsets. Once UMEIT is recast as a hy- 
perbolic problem, we face several difficulties. The equation is hyperbolic in the sense 
that one of the spatial variables plays the usual role of "time" in a second-order wave 
equation. Such a "time" variable has an orientation that depends on position x in X 
and also on the solution of the hyperbolic equation itself since the equation is nonlin- 
ear. Existence and uniqueness results for such equations need to be established, and we 
shall do so in sections [3] and [4] below adapting known results on linear and nonlinear 
hyperbolic equations that are summarized in [21, 29J. 

More damaging for the purpose of UMEIT and UMOT is the fact that hyperbolic 
equations propagate information in a stable fashion only when such information enters 
through a space-like surface, i.e., a surface that is more orthogonal than it is tangent to 
the direction of "time" . In two dimensions of space, the time-like and space-like variables 
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can be interchanged so that when n = 2, unwanted singularities can propagate inside 
the domain only through points with "null-like" normal vector, and in most settings, 
such points have (surface Lebesgue) zero-measure. In n = 2, it is therefore expected 
that spurious instabilities may propagate along a finite number of geodesies and that 
the reconstructions will be stable otherwise. 

In dimensions n > 3, however, a large part of the boundary dX will in general be 
purely "time-like" so that the information available on such a part of the surface cannot 
be used to solve the inverse problem in a stable manner [21J . Only on the domain 
of influence of the space-like part of the boundary do we expect to stably solve the 
nonlinear hyperbolic equation, and hence reconstruct the unknown conductivity cx(x). 

Special geometries and special boundary conditions. As we mentioned earlier, 
the partial reconstruction results described above can be improved in the setting of 
multiple measurements. Once several measurements, and hence several potential "time- 
like" directions are available, it becomes more likely that o can be reconstructed on 
the whole domain X. In the setting of well-chosen multiple measurements, the theories 
developed in [7J [H] indeed show that o can be uniquely and stably reconstructed on X. 

An alternative solution is to devise geometries of X and of the boundary conditions 
that guarantee that the "time-like" part of the boundary dX is empty. Information 
can then be propagated uniquely and stably throughout the domain. In section |3J we 
consider several such geometries. The first geometry consists of using an annulus-shaped 
domain and to ensure that the two connected components of the boundary are level-sets 
of the solution u. In such situations, the whole boundary dX turns out to be "space- 
like". Moreover, so long as u does not have any critical point, we can show that the 
reconstruction can be stably performed on the whole domain X. 

Unfortunately, only in dimension n = 2 can we be sure that u does not have any 
critical point independent of the unknown conductivity o. This is based on the fact that 
critical points in a two-dimensional elliptic equation are necessary isolated as used e.g. 
in pQ and our geometry simply prevents their existence. In three dimensions of space, 
however, critical points can arise. Such results are similar to those obtained in [12] in 
the context of homogenization theory and are consistent with the analysis of critical 
points in elliptic equations as in, e.g., JT3J [IS] . 

In dimension n > 3, we thus need to use another strategy to ensure that one vector 
field is always available for us to penetrate information inside the domain in a unique 
and stable manner. In this paper, such a result is obtained by means of boundary 
conditions / in (TP) that are "close" to traces of appropriate complex geometric optics 
(CGO) solutions, which can be constructed provided that a(x) is sufficiently smooth. 
The CGO solutions are used to obtain required qualitative properties of the solutions 
to linear elliptic equations as it was done in the setting of other hybrid medical imaging 
modalities in, e.g., 0, [9] dU [30]; see also the review paper [6]. 

The rest of the paper is structured as follows. Section [2] presents the derivation of 
the functional H(x) = cr(a;) |Vw| 2 from ultrasound modulation of a domain of interest 
and the transformation of the inverse problem as a nonlinear hyperbolic equation. In 
section [3j local results of uniqueness and stability are presented adapting results on 
linear hyperbolic equations summarized in [2H]. These results show that UMEIT and 
UMOT are indeed much more stable modalities than EIT and OT. The section concludes 
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with a local reconstruction algorithm, which shows that the nonlinear equation admits 
a solution even if the available data are slightly perturbed by, e.g., noise. The existence 
result is obtained after an appropriate change of variables from the result for time- 
dependent second-order nonlinear hyperbolic equations in [5TJ . Finally, in section |4j 
we present global uniqueness and stability results for UMEIT for specific geometries or 
specific boundary conditions constructed by means of CGO solutions. 



2 Derivation of a non-linear equation 

Ultrasound modulation. A methodology to combine high contrast with high reso- 
lution consists of perturbing the diffusion coefficient acoustically. Let an acoustic signal 
propagate throughout the domain. We assume here that the sound speed is constant 
and that the acoustic signal is a plane wave of the form pcos(k ■ x + <p) where p is the 
amplitude of the acoustic signal, k its wave-number, and ip an additional phase. The 
acoustic signal modifies the properties of the diffusion equation. We assume that such 
an effect is small but measurable and that the coefficient in ^ is modified as 

a E (x) — a(x)(l + ecos(k ■ x + <p)), (2) 

where e = pT is the product of the acoustic amplitude pel and a measure T > of the 
coupling between the acoustic signal and the modulations of the constitutive parameter 
m Q. For more information about similar derivations, we refer the reader to [2, 10, 23J. 

Let u be a solution of ([!]) with fixed boundary condition /. When the acoustic field 
is turned on, the coefficients are modified as described in (|2j) and we denote by u £ the 
corresponding solution. Note that U- £ is the solution obtained by changing the sign of 
p or equivalently by replacing ip by ip + tt. 

By the standard continuity of the solution to Q with respect to changes in the 
coefficients and regular perturbation arguments, we find that u £ = u + eu\ + 0(e 2 ). 
Let us multiply the equation for u e by w_ e and the equation for ti_ e by u £ , subtract the 
resulting equalities, and use standard integrations by parts. We obtain that 

f dii— du 
(<r s - o-_ £ )Vm £ ■ Vu_ £ dx = / o-_ e — —u £ - a £ —^u- £ da. (3) 
x Jax du dv 

Here, v{x) is the outward unit normal to X C IR n at x 6 dX and as usual ^ = v • V. 
We assume that a £ d v u £ is measured on dX, at least on the support of u E — f for all 
values e of interest. Note that the above equation still holds if the Dirichlet boundary 
conditions are replaced by Neumann (or more general Robin) boundary conditions. Let 
us define 

1 f du Oti 

J e := - / cr_ £ — -u £ - a £ — -U- £ da = eJ 1 + 0{e z ). (4) 

2 J dx dv dv 

The term of order 0(e 2 ) vanishes by symmetry. We assume that the real valued func- 
tions J\ = Ji(k,<p) are known. Such knowledge is based on the physical boundary 
measurement of the Cauchy data (u e , a £ d u u £ ) on dX. 

Equating like powers of e, we find at the leading order that: 



x 



[cr(x)Vuo • Vmq(x)] cos(A; ■ x + <p)dx = Ji(k, ip). (5) 
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This may be acquired for all k € M. n and <p = 0, |, and hence provides the Fourier 
transform of 

H{x) = cr(x)\Vu \ 2 (x). (6) 

Upon taking the inverse Fourier transform of the measurements (|5|, we thus obtain the 
internal functional ([6]). 



Nonlinear hyperbolic inverse problem. The forward problem consists of assuming 
a and f(x) known, solving ([T| to get u(x), and then constructing H{x) = cr(x)|Vu| 2 (x). 
The inverse problem consists of reconstructing a and u from knowledge of H(x) and 
fix). 

As we shall see, the linearization of the latter inverse problem may involve an op- 
erator that is not injective and so there is no guaranty that u and a can be uniquely 



reconstructed; see Remark 3.4 below. In this paper, we instead assume that the Neu- 
mann data ov ■ Vu and the conductivity <j(x) on OX are also known. We saw that 
measurements of Neumann data were necessary in the construction of H(x), and so our 
main new assumption is that a(x) is known on dX. This allows us to have access to 
v ■ Vu on dX. 

Combining ([I]) and ^ with the above hypotheses, we can eliminate a from the 
equations and obtain the following Cauchy problem for u(x): 

H ( x\ du 
- V ■ — A r Vu = in X, u = f and — = j on dX, (7) 
\\/u\ 2 (x) ov 

where (H,f,j) are now known while u is unknown. Thus the measurement operator 
maps (a, u) to (H,f,j) constructed from u(x) solution of ([T]). Although this problem 
(J7]) may look elliptic at first, it is in fact hyperbolic as we already mentioned and this 
is the reason why we augmented it with full (redundant) Cauchy data. In the sequel, 
we also consider other redundant measurements given by the acquisition of H(x) = 
a(x)\Vu\ 2 (x) for solutions u corresponding to several boundary conditions f(x). A 
general methodology to uniquely reconstruct u{x) from a sufficient number of redundant 
measurements has recently been analyzed in [71 [T4"] . 
The above equation may be transformed as 

(I - 2Vu ® V«) : V 2 m + V In if • Vw = in A, u = f and — = j on dX. (8) 



dv 



Here Vn = With 

|V«| 



gii = g^iVu) = -8 ij + 2(V«)»(Vw) j and k { = -(V\nH) h (9) 
then ^ is recast as: 

g lJ (Vu)d 2 i u + k l d i u = in A, u = f and — = j ondX. (10) 



dv 



Note that g %0 is a definite matrix of signature (l,n — 1) so that (10) is a quasilinear 
strictly hyperbolic equation. The Cauchy data / and j then need to be provided on a 
space-like hyper-surface in order for the hyperbolic problem to be well-posed [20] ■ This 
is the main difficulty when solving ^ with redundant Cauchy boundary conditions. 
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3 Local existence, uniqueness, and stability 



Once we recast ^ as the nonlinear hyperbolic equation ( |To| ) , we have a reasonable 
framework to perform local reconstructions. However, in general, we cannot hope to 
reconstruct u(x), and hence cr(x) on the whole domain X, at least not in a stable manner. 
The reason is that the direction of "time" in the second-order hyperbolic equation is 
Vu(x). The normal v{x) at the boundary dX separates the (good) part of dX that is 
"space-like" and the (bad) part of dX that is "time-like" ; see definitions below. Cauchy 
data on space-like surfaces such as t = provide stable information to solve standard 
wave equations where as in general it is known that arbitrary singularities can form in 
a wave equation from information on "time-like" surfaces such as x = or y = in a 
three dimensional setting (where (t,x,y) are local coordinates of X) |20j . 

The two-dimensional setting n = 2 is unique with this respect since the numbers of 
space-like and time-like variables both equal 1 and "t" and "x" play a symmetric role. 
Nonetheless, if there exist points at the boundary of dX such that v(x) is "light-like" 
(null), then singularities can form at such points and propagate inside the domain. As 
a consequence, even in two dimensions of space, instabilities are expected to occur in 
general. 

We present local results of uniqueness, stability of the reconstruction of u and a in 
These results are based on the linear theory of hyperbolic equations with 



section 3.1 



general Lorentzian metrics [29]. In section 3.2, we adapt results in [21] to propose a local 
theory of reconstruction of u(x), and hence a(x), by solving (10) with data (H, f,j) that 
are not necessarily in the range of the measurement operator (u, a) i— > (H, f,j), which 
to (u, a) satisfying ([I]) associates the Cauchy data (/, j) and the internal functional H. 



3.1 Uniqueness and stability 

Stability estimates may be obtained as follows. Let (u, a) and (u, a) be two solutions of 
Q and the Cauchy problem (10) with measurements (H,f,j) and (H,f,j). Note that 
after solving (10), we then reconstruct the conductivities with 



H 



H 



(7lX 



(T(X) 



Vu 



r,\2 



x 



The objective of stability estimates is to show that (u — u, o — a) are controlled by 
(H — H, f — f,j — j), i.e., to show that small errors in measurements (that are in the 
range of the measurement operator) correspond to small errors in the coefficients that 
generated such measurements. 

Some algebra shows that v = u — u solves the following linear equation 



H 



ii 1 2 



Vu <g> {Vu + Vu) 
IVwP 



Vv 



H -H. 



Vu 



r, 12 



-Vu 



0. 



with Cauchy data / — / and j — j, respectively. Changing the roles of u and u and 
summing the two equalities, we get 



V-l 



H 



\Vu\ 2 \Vu\ 2 



(Vu+Vu)®(Vu+Vu)-(\Vu\ 2 +\Vu\ 2 )l}Vv+5H 



Vu Vu " 
|Wl 2 + |Vd 2 . 



0. 
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The above operator is elliptic when Vu • Vu < and is hyperbolic when Vu • Vu > 0. 
Note that Vu • Vu > on dX when j — j and f — f are sufficiently small. We obtain a 
linear equation for v with a source term proportional to 5H = H — H . For large amounts 
of noise, Vu may significantly depart from Vu, in which case the above equation may 
loose its hyperbolic character. However, stability estimates are useful when 6H is small, 
which should imply that u and u are sufficiently close, in which case the above operator 
is hyperbolic. We assume here that the solutions u and u are sufficiently close so that 
the above equation is hyperbolic throughout the domain. We recast the above equation 
as the linear equation 

S v (x)^ j v + VdiU + d i (pSH) = mX, v = f-f, ^ = j - j on dX, (12) 

for appropriate coefficients Q 13 , F and l % . Now q %3 is strictly hyperbolic in X (of signature 
(1, n — 1)) and is given explicitly by 



Q(x) = ,J (Vu + Vu) ® (Vu + Vu) - (\Vu\ 2 + \Vu\ 2 )l\ 

\yu\ z \yu\ z I J 

= a(x) ( e(x) ® e(x) — fi 2 (x) (I — e(x) ® e(x)) 



(13) 



where 

, Vu + Vu 2 |V«| 2 + |V£t| 2 

e x = . J, /? 2 x = ' ' ' ' x , 14 

|Vw + Viz | | V« + \7u\ z — (\ Viip + | VitpJ 

and a(x) = ^~^ u ^ (\Vu + Vu\ 2 — (|Vu| 2 + |Vtt| 2 )) is the appropriate (scalar) normal- 
ization constant. Here, e(x) is a normal vector that gives the direction of "time" and 
/3(x) should be seen as a speed of propagation (close to 1 when u and u are close). When 
e is constant, then the above metric, up to normalization, corresponds to the operator 

We also define the Lorentzian metric f) = so that Ijy are the coordinates of the 
inverse of the matrix q 13 . We denote by (-, •) the bilinear product associated to f) so that 
(u,v) = \)ijU % v 3 where the two vectors u and v have coordinates u % and v\ respectively. 
We verify that 

f)(x) = ^y(e(x)®e(x)-^y(/-e(x)®e(x))). (15) 



The main difficulty in obtaining a solution v to (12) arises because v{x) is not time- 
like for all points of dX. The space-like part S 9 of dX is given by the points x G dX 
such that v(x) is time-like, i.e., f)(z/(x), u(x)) > 0, or equivalently 

\u(x)-e(x)\ 2 > 1 + ^ 2(x) xedX. (16) 

Above, the "dot" product is with respect to the standard Euclidean metric and v is a 
unit vector for the Euclidean metric, not for the metric f). The time-like part of dX is 
given by the points x G dX such that l)(z/(x), ^(x)) < (i.e., z/(x) is a space-like vector) 
while the light-like (null) part of dX corresponds to x such that l)(z/(x), v{x)) = (i.e., 
v{x) is a null vector). 
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When j = j on dX so that Vm(i) = Vtt(x) and (3(x) = 1 for x G <9X (see also the 
proof of Theorem 3.1| below), then the above constraint becomes: 



\u{x) ■ Vu(x)\ 2 > - xedX. (17) 

In other words, when such a constraint is satisfied, the differential operator is strictly 
hyperbolic with respect to u(x) on E 5 . Once E 9 is constructed, we need to define 
its domain of influence X g C X, i.e., the domain in which v can be calculated from 
knowledge of its Cauchy data on E g . In order to do so, we apply the energy estimate 
method for hyperbolic equations described in [29j Section 2.8]. We need to introduce 
some notation as in the latter reference; see Figure [TJ 




Figure 1: Construction of the domain of influence O (hatched area in the above picture). 
The unit vectors e indicate the "time" direction of the Lorentzian metric f). The surface 
S 2 (s) has a normal vector z/ 2 (x) that forms a sufficiently small angle with e so that E 2 (s) 
is a space-like surface, as is Ex C E 9 with an angle such that \v\ ■ e| is also sufficiently 
close to 1. 

Let Ex be an open connected component of E fl . We assume here that all coefficients 
and geometrical quantities are smooth. By assumption, Ei is space-like, which means 



that the normal vector v\ is time-like and hence satisfies (16). Let now S 2 (s) C X be 



a family of (open) hyper-surfaces that are also space-like with unit (with respect to the 



Euclidean metric) vector i^ 2 (x) that is thus time-like, i.e., verifies (16). We assume that 



the boundary of S 2 (s) is a co-dimension 1 manifold of Ex- Let then 

0(s) = |J E 2 (r), 



0<T<S 



which we assume is an open subset of X. In other words, we look at domains of influence 
0(s) of Ei that are foliated (swept out) by the space-like surfaces S 2 (r). Then we have 
the following result: 
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Theorem 3.1 (Local Uniqueness and Stability.) Let u and u be two solutions of 
(|7| sufficiently close in the W 1 '°°(X) norm and such that \Vu\, \Vu\, H and H are 
bounded above and below by positive constants. This ensures that Q constructed in (13) 



is strictly hyperbolic and that a(x) and (3(x) in (14) are bounded above and below by 
positive constants. 

Let Si be an open connected component of S 9 the space-like component of dX and 
let the domain of influence O = O(s) for some s > be constructed as above. Let us 
define the energy 

E{dv) = (dv,u 2 ) 2 - -(dv,dv)(v 2 ,v 2 ). (19) 

Here, dv is the gradient of v in the metric f) and is thus given in coordinates by Q^djV. 
Then we have the local stability result 



E(dv)dx<C( / \f-f\' 

o V JSi 



\j-]\ 2 da + / \V5H\ 2 dx 



o 



(20) 



where dx and da are the standard (Euclidean) volume and (hyper) surface measures on 
O and Si, respectively. 

The above estimate is the natural estimate for the Lorentzian metric h. For the 
Euclidean metric, the above estimate may be modified as follows. Let u 2 (x) be the unit 
(for the Euclidean metric) vector to x G E 2 (s) and let us define c(x) := v 2 (x) ■ e(x) with 



e(x) as in (14). Let us define 



6 :-- 



mm 



c 2 (x) 



1 



l + /3 2 (x) 



(21) 



We need 9 > for the metric f) to be hyperbolic with respect to u 2 (x) for all x G O. 
Then we have that 



\v 2 \ + \Vv\ 2 + (d - a) 1 dx < 



o 



C 



\f-f\ 2 + \j-j\ 2 da+ \VSH\ 2 dx), (22) 



o 



where a and a are the reconstructed conductivities given in (11). Provided that data 
are equal in the sense that f = f , j = j , and H = H , we obtain that v = and the 
uniqueness result u = u and o = a. 

Proof. That I) is a hyperbolic metric is obtained for instance if u and u are sufficiently 
close in the W^°°{X) norm and if |Vw|, |Vw|, H and H are bounded above and below 



by positive constants. The derivation of (20) then follows from [29j Proposition 8.1] 
using the notation introduced earlier in this section. The volume and surface measures 
dx and da are here the Euclidean measures and are of the same order as the volume 



and surface measures of the Lorentzian metric h. This can be seen in (15) since a and 
(3 are bounded above and below by positive constants. 



Then (20) reflects the fact that the energy measured by the metric f) is controlled. 



However, such an "energy" fails to remain definite for null-like vectors (vectors v such 
that h(v,v) = 0) and as x approaches the boundary of the domain of influence of S 9 , 
we expect the estimate to deteriorate. 

Let x G O be fixed and define v = v 2 {x) and e = e(x). Let us decompose v = 
ce + s'e- 1 , where ce is the orthogonal projection of v onto e and s'e 1 - := v — ce the 
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projection onto the orthogonal subspace of IR ra with e 1 - a unit vector. For a vector 
v = fie + v^e 1 - + w' (standing for dv) with w' orthogonal to e and e 1 - (and thus 
vanishing if n = 2), we need to estimate 



1 1 



E(v) = tf(v. u) - ^(v, v)h(u, v) = - ( Vl c - v 2 sf -\{vl- (v 2 2 + |w| 2 ))(c 2 - s 2 ) 

where we have conveniently defined v% = P~ 1 v' 2 , w = /3 _1 w', and s = After some 

straightforward algebra, we find that 

£( v ) = 4Jr^l w l 2 + - 2 (kc 2 + s 2 )(v 2 1 +vl)-2v 1 v 2 cs) 
a z p z or V2 J 

Since /3 is bounded above and below by positive constants, we need to bound (c — s) 
from below or equivalently (/3c — s') 2 from below. Some algebra shows that 

a / 2 1 / 3c + s Vfl ^ 

^ C -TT^ = TT^ (/3c - s) - 



Since 6* < 1, this shows that 



EM > C9 2 W 



2 1 |2 



for a constant C that depends on the lower and upper bounds for and a but not on 
the geometry of v. Note that the behavior of the energy in 9 2 is sharp as the bound is 
attained for v\ = v% with w = 0. This proves the error estimate for v = Vf . Since v is 
controlled on Ei, we obtain control of v on O by the Poincare inequality. Now a — a is 
estimated by H — H and by V« — V£t = Vf and hence the result. 

In other words, the angle (j>(x) between e(x) and u 2 (x) must be such that /3(x) — 
tan0(x) > 9 2 in order to obtain a stable reconstruction. When 5H is small, then 
Vu — Vu is small so that is close to 1. As a consequence, we obtain that the constraint 
of hyperbolicity of f) is to first order that tan0(x) < 1, which is indeed the constraint 



(16) that holds when Vw = Vm on dX. 

For the uniqueness result, assume that u and u are two solutions of ([7]). We define 
e(x) = Vm and /3 2 = 1. Then v = on Si implies by the preceding results that v — 
in a vicinity of Si in O so that u = u in the vicinity of Si. This shows that u = u in O 
and hence in all the domain of dependence of S 9 constructed as above. □ 

Let us conclude this section with a few remarks. 



Remark 3.2 In two dimensions, we can exchange the role of space-like and time-like 
variables since both are one-dimensional and find, at least for sufficiently simple ge- 
ometries, that the complement of the domain of influence of S 9 in X is the domain of 
influence of the complement of S g in OX. We thus obtain stability of the reconstruction 
in all of X except in the vicinity of the geodesies for the metric g that emanate from 
dX in a direction v(x) that is null-like, i.e., a vector such that i)(y(x),u(x)) = 0, or 
equivalently such that \u{x) ■ e(x)| 2 = |. 

In three (or higher) dimensions, however, no such exchange of the role of time and 
space is possible. All we can hope for is a uniqueness and stability result in the domain 
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Figure 2: Geometry of the domain of influence in Euclidean geometry with a = 1 and 
u = X\ on a domain X given by an ovoid. In a three dimensional geometry, we can 
consider the above picture a cross section at y — of a three dimensional domain of 
revolution about the axis ei = Vw. The vector u(x) is a "null-vector" making an angle 
at 45 degrees with Vm. In two dimensions, S 9 is the union of two connected components 
where as in thee dimensions, £ 5 is composed of a unique connected component in dX. 
The hatched domain corresponds to X\X g , the part of the domain X that is not the 
domain of influence of S g . In two dimensions, Vu -1 = e 2 may also play the role of "time" 
so that X\X g is the domain of influence of E\E ff . In three dimension, the hatched region 
is not accessible with the techniques developed in this paper. 



of influence of £ 9 . The solution v and the conductivity a are not stably reconstructed on 
the rest of the domain without additional information, coming e.g. from other boundary 
conditions f(x). The case of redundant measurements of this type is considered in 
section 



4.2| below, and is analyzed in a different context in [TJ [H] . 



Remark 3.3 Assuming that the errors on the Cauchy data / and j are negligible, we 
obtain the following stability estimate for the conductivity 

C 

\W - v\\lho) < - H\\ H i( X ). (23) 

The measurements are of the form H(x) = cx(x)\Vu\ 2 (x), which imposes reasonably 
restrictive assumptions on a ensuring that Vm is a solution in H 2 (Q). Under additional 
regularity assumptions on cr, for instance assuming that H e H S (X) for s > 2, we find 
that 

\W - *\\ LH o) <j\\H- Hlf-fa \\H + H\\l s{x) , (24) 

by standard interpolation. We thus obtain a standard Holder estimate in the setting 
where the error in the measurements is quantified in the square integrable sense. 
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Remark 3.4 The linearization of ^ in the vicinity of o"o = 1 with only Dirichlet data 
is an ill-posed problem when X is a two-dimensional disc. Indeed, assume Dirichlet data 
of the form f(x) = X\ in Q so that the unperturbed solution is u = x% in X. This 



shows that e(x) = e 1 in the definition (13) so that () = ei <8> ei — e 2 <8> e 2 in (15). In 



other words, the linearized problem consists of solving 
d 2 u d 2 u 



dx\ 2 dx2 2 



in X = {xf + x 2 2 < 1}, u = f on dX. 



The general solution to the above equation is of the form F(x\ — x 2 ) + G(x\ + x 2 ) and 
there is an infinite number of linearly independent solutions to the above equation with 
/ = 0. The linearization of the UMEIT problem in this specific geometry provides an 
operator that is not injective. 

3.2 Reconstruction of the conductivity 

The construction of the solution u, from which we deduce the reconstruction of a(x), 



requires that we solve the nonlinear equation (10). Let us assume that g tJ is given as in 
([9]) and that the vector field h and the source terms / and j are smooth given functions. 
Then we can construct a unique solution to (10) locally in the vicinity of the part of 
dX that is space-like. In this section, we assume that the geometry and the coefficients 
of the wave equation are sufficiently smooth. 

Let a point x° G E s , the space-like part of dX, i.e., so that g(v(x°), z/(x )) > 77 > 
0. In the vicinity of x°, which we now call 0, we parameterize dX by the variables 
(yi, . . . ,y n -i) an d denote by yo the signed distance to dX. In the vicinity of x° = 0, 
y — F(x) is a diffeomorphism from a neighborhood U of x = to the neighborhood 
V = F(U) of y = 0. Moreover, locally, DF is close to the identity matrix (after 
appropriate rotation of the domain if necessary) if U is sufficiently small. We denote by 
Jp = det(DF) the Jacobian of the transformation. 

Let us come back to the equation 

- V ■ a(x)Vu = -V ■ — J , Vu = in X u = f and — = j on dX. (25) 

\\/u\ 2 (x) av 

We define v(y) = u(x), i.e., v = F*u and then verify that (Vu)(x) = DF t oF' 1 (y)'Vv(y). 
In the y coordinates, we find that 

-V-F*aVv = 0, F(U), 

where we have the standard expression in the y coordinates: 

F.a(y) = a(y)DFDF t o F~\y), a = J F l a o F~ l . 
We may thus recast the above equation as the following non-linear equation 

~ DFDF 1 o F^ 1 

- V ■ H lDFtVv{2 = 0, H — F*{J F l H) = J?H o F . (26) 

Note that the boundary conditions are now posed on the surface yo = where 
v(0,y') = FJ(0,y') and d yo v(0,y') = a{y')FJ(0,y') + 0(y')F*j(O,y'), 
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with a close to and (3 close to 1 on V = F{U). It remains to differentiate in (26) to 
obtain after straightforward but tedious calculations the expression 

gp 1 df j v + hpd i v = 0, F(U), (27) 

with the above "initial" conditions at yo = 0, where 

g% = -(DFDF%V k + 2(DFDF^v)\DF*Vvy i 

h F = -(VlnF • DFDF*)* - (V ■ DFDF 1 ) 1 + 2{DFDF^v) i {DF^Vv)d^DFi. 

_J28) 

When F = I, we recover ([8]). The non-linear terms now involve functions of DF tJ \7v. 

Note that gp = DF g DF l if we denote DF^v = V« and thus transforms as a tensor 
of type (2,0). As a consequence, the metric (a tensor of type (0,2)) gp 1 = DF~ 1 gDF~ t 
since g' 1 = g as can be easily verified. Let Vp = F*v = DFv o F" 1 be the push-forward 
of the normal vector seen as a vector field. At x°, the change of variables is such that 

d d 

9KF(xO)(gy-^ Qy-) = 9p,p (x0) (VF, V F ) = (v, v) = g x o{u, V)>T)>0. 

This shows that gp remains hyperbolic in the vicinity of F(x°) since DFDF 1 = I at 
y = by construction and DF is smooth. Moreover, the above is equivalent to 

£)2 n—1 
3= 



dy°dy° y ' y ^dyidyi 



where IT 15 ] — T 2- ^ ne a b° ve is nothing but the fact that gp is hyperbolic in the 

vicinity of y — 0. Note that 7 1 - 7 = 7 u '(a;, DF l Vv). 

We thus have a nonlinear hyperbolic equation of the form 

(□ + j ij (x, DF^v)dij + h^x, DF^Vv)di)v = 0, y > 0, y' e M n_1 , , 

Since propagation in a wave equation is local we can extend the boundary conditions 
for y = (0,y') outside the domain F(U) by v o = and d yo v = 1 and the functions 7 U 
and h l by outside of F{U). This allows us to obtain an equation posed on the half 
space yo > 0. 

The nonlinear functions j^^x, DF^v) and h l (x, DF^v) are smooth functions of 
Vf except at the points where Vt> = 0. However, we are interested in solutions such 
as Vf does not reach to preserve the hyperbolic structure of g^ . Note that |Vf | is 
bounded from below by a positive constant on y = by assumption. We obtain a bound 
on the uniform norm of the Hessian of v, which implies that at least for a sufficiently 
small interval y G (0, t ), \Vv\ does not vanish and ^ and h % can then be considered 
as smooth functions of x and Vv . 

Using Theorem 6.4.11 in [21] and the remark following (6.4.24) of that reference, the 
above equation satisfies the hypotheses to obtain an a priori estimate for 

M (Vo)= ^ 11^(2/0,0 lU^Rn-l), 
|a|<re+2 
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with k the smallest integer strictly greater than By Sobolev imbedding, this implies 
that the second derivatives of v are uniformly bounded so that for at least a small 
interval, |Vt>| is bounded away from 0. 

Once v, and hence u is reconstructed, at least in the vicinity of the part E g of OX 
that is space like for Vw, we deduce that 

, v H(x) 



\Vu\ 2 (x) 



Note that V« cannot vanish by construction so that the above equality for cr(x) is well- 
defined. We already know that a solution to the above nonlinear equation exists in the 
absence of noise since we have constructed it by solving the original linear equation. In 
the presence of significant noise, the nonlinear equation may behave in a quite different 
manner than that for the exact solution. However, the above construction shows that 
the nonlinear equation can be solved locally if the measurement H(x) is perturbed by 
a small amount of noise. 



4 Global reconstructions of the diffusion coefficient 

The picture in Fig. [2] shows that in general, we cannot hope to obtain a global recon- 
struction from a single measurement of H(x) even augmented with full Cauchy data. 
Only Cauchy data on the space-like part of the boundary can be used to obtain stable 
reconstructions. 

Global reconstructions have been obtained from redundant measurements of the form 
Hij = Si ■ Sj with Si = ^/aVui and Ui solution of (fij) with Dirichlet conditions / = /j, in 
[13] in the two dimensional setting and in [7] in the two- and three- dimensional settings; 
see also [23] . 

This section analyzes geometries in which a unique measurement H(x) or a small 
number of measurements of the form H(x), augmented with Cauchy data (/, j) allow one 
to uniquely and stably reconstruct a(x) on the whole domain X. These reconstructions 
are obtained by (possibly) modifying the geometry of the problem so that the domain 
where cr(x) is not known lies within the domain of dependence of E s . We consider 
two scenarios. In the first scenario, considered in section |4.1[ we slightly modify the 
problem to obtain a model with an internal source of radiation /. Such geometries 
are guaranteed to provide a unique global reconstruction in dimension n = 2 but not 
necessarily in higher spatial dimensions, where global reconstructions hold only for a 
certain class of coefficients a(x). In the second scenario, analyzed in section 4.2, we 



consider a setting where reconstructions are possible when the Lorentzian metric is the 



Euclidean (Lorentzian) metric, i.e., a = ft = 1 in (15). We then show the existence of 
an open set of illuminations / for three different measurements of the form H(x) such 
that the global result obtained for the Euclidean metric remains valid for arbitrary, 
sufficiently smooth coefficients a(x). 

4.1 Geometries with an internal source 

From the geometric point of view, the Cauchy data are sufficient to allow for full recon- 
structions when S 9 = dX, so that the whole boundary dX is space-like for the metric g, 
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and X is the domain of dependence of S g . This can happen for instance when dX is a 
level set of u and the normal derivative of u either points inwards or outwards at every 
point of dX. When X is a simply connected domain, the maximum principle prevents 
one from having such a geometry. However, when X is not simply connected, then such 
a configuration can arise. We will show that such a configuration (with X the domain 
of dependence of S g ) is always possible in two dimensions of space. When n > 3, such 
configurations hold only for a restricted class of conductivities <j(x) for which no critical 
points of u(x) exist. 

Let us consider the two dimensional case n — 2. We assume that X is an open 
smooth domain diffeomorphic to an annulus and with boundary dX = dX U dXi, see 
Fig. |3j We assume that / = on the external boundary OXq and / = 1 on the internal 
boundary dX\. The boundary of X is composed of two smooth connected components 
that are different level sets of the solution u to (JT1), which is uniquely defined in X. 




Figure 3: Geometry of an annulus in two space dimensions with boundaries S and S 1; 
the level sets where u = and u — 1, respectively^ The four curves correspond to four 
integrals of the flow of the gradient vector field Vw. 

In practice, such a domain X may be constructed as follows. As we do in the 
geometry depicted in Fig. [5] below, we embed X, the domain where a is unknown, into 
a larger domain X with, e.g., a(x) = <jq on X\X and with a hole where we impose the 
aforementioned boundary conditions. Then we have the following result: 

Proposition 4.1 Let X be the geometry described above with n = 2 and u(x) the 
solution to ([T]). We assume here that both the geometry and cr(x) are sufficiently smooth. 
Then \Vu\ is bounded from above and below by positive constants. The level sets E c = 
{x G X, u(x) = c} for < c < 1 are smooth curves that separate X into two disjoint 
subdomains. 

Proof. The proof of the first part is based on the fact that critical points of solutions 
to elliptic equations in two dimensions are isolated pQ. First of all, the Hopf lemma [T3] 
ensures that no critical point exists on the smooth closed curves E and Si. Let Xj be 
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the finite number of points where Vti(xj) = 0. At each Xi, the level set of u with value 
< Ci = u(xi) < 1 is locally represented by n« (n^ even) smooth simple arcs emanating 
from Xi that make an angle equal to 2ir/ni at x^ pp. For instance, if only two simple arcs 
emanate from xq, then these two arcs form a continuously differentiable curve in the 
vicinity of x$. Between critical points, level sets of u are smooth by the inverse function 
theorem. 

Let us assume that there is a point Xi with more than two simple arcs leaving x^. 
Let 7j, 1 < j < 4, be such arcs. If 71 meets another critical point, we pick one of the 
possible other arcs emanating from this critical point to continue the curve 71. This is 
always possible as critical points always have an even number of leaving simple arcs. 
The curve 71 cannot meet So or Si and therefore must come back to the point Xj. Let 
us assume the existence of a closed sub-loop of 71 that does not self-intersect and does 
not wind around Si (i.e., is homotopic to a point). In the interior of that close sub-loop, 
u is then constant by the maximum principle and hence constant on X by the unique 
continuation theorem [19] . This is impossible and therefore 71 must wind around S x . 
Let us pick a subset of 71, which we still call 71 that winds around Si once. The loop 
meets one of the other jj to come back to Xi, which we call 72 if it is not 71. Now let 
us follow 73. Such a curve also has to come back to Xi. By the maximum principle 
and the unique continuation theorem, it cannot come back with a sub-loop homotopic 
to a point. So it must come back also winding around S x . But 71 and 73 are then 
two different curves winding around Si. This implies the existence of a connected (not 
necessarily simply connected) domain whose boundary is included in 71 U 73. Again, 
by the maximum principle and the unique continuation theorem, such a domain cannot 
exist. So any critical point cannot have more than two simple arcs of level curves of u 
leaving it. 

So far, we have proved that any critical point exactly two arcs leaving Xj at an 

angle equal to 7r since by the maximum principle, critical points can not be local mimina 
or maxima. These two arcs again have to meet winding around Si. This generates what 
we call a single curve 71 with no possible self-intersection. Moreover, since all angles at 
critical points are equal to tt, the curve 71 is of class C 1 and piecewise of class C 2 . Let 
X c be the annulus with boundary equal to S1U71. On X c , u satisfies an elliptic equation 
with values u = 1 on Si and < u = c$ < 1 on 71. Since 71 is sufficiently smooth now 
(smooth on each arc with matching derivatives on each side of each critical point), it 
satisfies the interior sphere condition and we can apply the Hopf lemma [T7J Lemma 
3.4] to deduce that the normal derivative of u on 71 cannot vanish at Xi or anywhere 
along 71. There are therefore no critical points of u in X. By continuity, this means 
that |Vu| is uniformly bounded from below by a positive constant. Standard regularity 
results show that it is also bounded from above. 

Now let < c < 1 and S c be the level set where u = c. Such a level set separates 
X into two subdomains where < u < c and c < u < 1, respectively, by the maximum 
principle. We therefore obtain a foliation of X into the union of the smooth curves S c 
for < c < 1. Now let x G S c and consider the flow of Vw in both directions emanating 
from x. Then both curves are smooth and need to reach the boundary at a unique 
point. Since any point on So is also mapped to a point on Si by the same flow, this 
shows that S c is diffeomorphic to S and Si. □ 

The result extends to higher dimensions provided that |Vw| does not vanish with 
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exactly the same proof. Only the proof of the absence of critical points of u was purely 
two-dimensional. In the absence of critical points, we thus obtain that e(x) = Vu = u(x) 



so that u(x) is clearly a time-like vector. Then the local results of Theorem 3.1 become 
global results, which yields the following proposition: 

Proposition 4.2 Let X be the geometry described above in dimension n > 2 and u(x) 
the solution to ([I]). We assume here that both the geometry and <r{x) are sufficiently 
smooth. We also assume that \ Vu\ is bounded from above and below by positive constants. 



Then the nonlinear equation (10) admits a unique solution and the reconstruction of u 



and of a is stable in X in the sense described in Theorem 3.1 



Remark 4.3 The above geometry with a hole is not entirely necessary in practice. 
Formally, we can assume that the hole with boundary Si shrinks and converges to 
a point xq G dX at the boundary of the domain. Thus, the illumination / is an 
approximation of a delta function at xq. The level sets of the solution are qualitatively 
similar to the level sets in the annulus. Away from xq, the surface dX is a level set of the 
solution u and hence the normal to the level set is a time-like vector for the Lorentzian 
metric with direction e(x) = u(x). Away from xq, we can solve the wave equation 
inwards and obtain stable reconstructions in all of X but a small neighborhood of xq. 
This construction should also provide stable reconstructions in arbitrary dimensions 
provided that u does not have any critical point. 

In dimensions n > 3, however, we cannot guaranty that u does not have any critical 
point independent of the conductivity. If the conductivity is close to a constant where 
we know that no critical point exists, then by continuity of u with respect to small 
changes in a(x), u does not have any critical point and the above result applies. In 
the general case, however, we cannot guaranty that Vw does not vanish and in fact can 
produce a counter-example using the geometry introduced in [12] (see also [25J for the 
existence of critical points of elliptic solutions): 

Proposition 4.4 There is an example of a smooth conductivity such that u admits 
critical points. 

Proof. Consider the geometry in three dimensions depicted in Fig. |4j The domain X 
is a smooth, convex, domain invariant by rotation leaving e z invariant and by symmetry 
z — > —z and including two disjoint, interlocked, tori T\ and T 2 . The first torus T x is 
centered at c\ = (0, 0, 1) with base circle {e z + 2e x + a(cos4>e x + sin^e^), < 4> < 2ti} 
rotating around c\ in the plane (e x ,e z ) (top torus in Fig. 0) for a = |, say. The 
second torus T 2 is centered at c 2 = (0, 0, —1) with base circle {— e 2 + 2e y + a(cos0e x + 
sm(fre y ), < 4> < 2n} rotating around c 2 in the plane (e^, e z ) (bottom torus in Fig. [1]). 
We consider the boundary condition u = z on dX. 

We assume that a(x) = 1 + \(p(x) in ([!]), where f(x) is a smooth, non-trivial, non- 
negative function with non-vanishing support inside each of the tori 7\ and T 2 that 
respects the invariance by rotation and the symmetries of the two tori. We normalize 
ip(x) by 1 on the circles {e z + 2(cos (pe z + sin^e^), < 4> < 2tt} and {— e 2 + 2(cos <fie z + 
sin^e^), < <f> < 2tc} at the center of the volumes delimited by the two tori. When 
A = so that a(x) = 1, then u = z is the solution of the problem (jl]). As A, an hence a 
inside the tori, converges to +oo, the solution u is such that u converges to a constant 
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Figure 4: Geometry of a critical point. X is the ball of radius 4. The interlocked tori 
are the top torus 7\ and the bottom torus T 2 . 

C\ > on the support of cp inside 7\ and C 2 < on the support of <p inside T 2 . For 
A sufficiently large, by continuity of the solution u with respect to a, we obtain that 
w(0,0, 1) < since (0,0,1) is inside T 2 and w(0,0, — 1) > since (0,0,-1) is inside 
T\. Since the geometry is invariant by symmetry x — > —x and y — > —y, then so is 
the solution x and hence d x u(0,0, z) = d y u(0,0,z) = for all (0, 0, z) G X. Now the 
function z — > u(0, 0, z) goes from negative to positive to negative back to positive values 
as z increases, and so has at least two critical points. At these points, Vu = and 
hence the possible presence of critical points in elliptic equations in dimensions three 
and higher. □ 

Note that the above symmetries are not necessary to obtain critical points, which 
appear generically in structures of the form of two interlocked rings with high conduc- 
tivities as indicated above. We describe this results at a very intuitive, informal, level. 
Indeed, small perturbations of the above geometry and boundary conditions make that 
the level sets E c = {u = c} for c sufficiently large and c sufficiently small are simply 
connected co-dimension 1 manifolds with boundary on dX. When o is sufficiently large, 
then u converges to two different values c\ and c 2 inside the two discs (say one positive 
in T\ and one negative in T 2 ). Thus for o sufficiently large, the level set u — Ci, assuming 
it does not have any critical point, is a smooth locally co-dimension 1 manifold, by the 
implicit function theorem, that can no longer be simply connected. Thus as the level 
sets c decrease from high values to ci, they go through a change of topology that can 
only occur at a critical point of u [26J. 
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4.2 Complex geometric optics solutions and global stability 

Let us now consider a domain X where cr(x) is unknown and close to a constant a . Let 
us assume that X is embedded into a larger domain X and that we can assume that 
cr(x) is known and also close to the constant cr . Then, it is not difficult to construct 
X so that X lies entirely within the domain of dependence of £ g ; see for instance the 
geometry depicted in Fig. |5j 




Figure 5: Geometry of a domain where the reconstruction of the unknown cr on X is 
possible from a single measurement. The geometry of the Lorentzian metric is repre- 
sented when cr(x) = cq. By continuity, the domain of influence of £ 9 includes X for all 
smooth conductivities a(x) sufficiently close to a . 

For the rest of the section, we show that global reconstructions can be obtained 
for general sufficiently smooth metrics provided that three well-chosen measurements 
are available. This result is independent of spatial dimension. The measurements are 
constructed by means of complex geometrical optics solutions. 

Let k be a vector in W 1 and k L be a vector orthogonal to k of same length. Let 
p = ik + k L be a complex valued vector so that p ■ p = 0. Thus, e p ' x is harmonic and 
Ve p ' x = pe p ' x . The latter gradient has a privileged direction of propagation p, which is, 
however, complex valued. Its real and imaginary parts are such that 

e -k± X Q Ve p-x = \k\8(x), e- k± - x me p - x = \k\6^{x), (30) 

where 9{x) = k cos k-x + k 1 - sin k-x and 9 ± (x) = —k sin k-x + k 1 - cos k-x. As usual, k — pr. 

Consider propagation with Cauchy data given on a hyperplane with normal vector 
k G § n ~ 1 . We want to make sure that we always have at our disposal a Lorentzian metric 
for which A; is a time-like vector so that the available Cauchy data live on a space-like 
surface for that metric. For the rest of the section, we assume that k = \k\e\ and that 
k 1 - = \k\e 2 so that 

9{x) = k cos \k\xi + k 1 sin \k\xi and 9 ± (x) = — k sin \k\x\ + k 1 cos \k\x\ (31) 

For a vector field with unit vector 9(x), we associate the Lorentz metric with direction 
9 given by t) e = 29 <g) 9 - I. 
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The Lorentzian metrics with directions 8(x) and 9 ± (x) oscillate with x±. A given 
vector k therefore cannot be time-like for all points x. However, we can always construct 
two different linear combinations of these two directions that form time-like vectors for 
a given range of k ■ x = \k\x\. Such combinations allow us to solve the wave equation 
forward and obtain unique and stable reconstructions on the whole domain X. The 
above construction with e p ' x harmonic can be applied when cr(x) = <Jq a constant. It 
turns out that we can construct complex geometric optics solutions for arbitrary, suffi- 
ciently smooth conductivities a(x) and obtain global existence and uniqueness results 
in that setting. We state the following result. 

Theorem 4.5 Let o be extended by = 1 on IR n \X ; where X is the domain where a 
is not known. We assume that a is smooth on W 1 . Let a(x) — l be supported without loss 
of generality on the cube (0, 1) x (— |, \) n ~ 1 ■ Define the domain X = (0, 1) x B n ^i(a), 
where B n -x(a) is the n — 1-dimensional ball of radius a centered at and where a is 
sufficiently large that the light cone for the Euclidean metric emerging from B n _i(a) 
strictly includes X. Then there is an open set of illuminations (fi,f 2 ) such that if u\ 
and U2 are the corresponding solutions of ([!]), then the following measurements 

H 11 (x) = a(x)\Vu 1 \ 2 (x), H 22 (x) = a{x)\Vu 2 \ 2 {x), H l2 {x) = a(x)\V( Ul + u 2 )\ 2 (32) 

with the corresponding Cauchy data (/i,ji), (f 2 ,j 2 ) and (fx + f 2 ,ji + 32) at Xi = 
uniquely determine a(x). Moreover, let be measurements corresponding to a and 
Ji) an d (/2> 32) the corresponding Cauchy data at x\ = 0. We assume that a(x) — 1 
and a{x) — 1 (also supported in (0, 1) x (— |, \) n ^ 1 ) are smooth and such that their norm 
in H 2+ 3+e (]R n ) for some e > are bounded by M . Then for a constant C that depends 
on M, we have the global stability result 

\W~ °\\l*(x) < Cy\dc- dcWiL^Bn-^aW + WVHij-VHijWvux))- (33) 

Here, we have defined I = {(1, 1), (1, 2), (2, 2)} and dc = (/1, ji, /2, 32) with dc being 
defined similarly. 

Proof. We recall that k = \k\e\ and k x = \k\e 2 . The proof is performed iteratively 
on layers ti-\ < x\ < U with ti — jr for < i < N and iV = N(k) (to be determined) 
sufficiently large but finite for any given sufficiently smooth conductivity o~(x). Here, 
k = \k\ei is the vector in IR n used for the constructions of the CGO solutions. We define 
yi — (ti, 0, . . . , 0) for < i < N. Define two vectors close to e x as 

p = wei + vT— w 2 e 2 , q = we\ — \/T— w 2 e 2 , 

with w < 1 sufficiently close to 1 such that the light cones (for the Euclidean metric) 
emerging from B n _i(a) for the Lorentzian metric with main directions p and q still 
strictly include X; see Fig. [6] All we need is that the radius a be chosen sufficiently 
large so that any Lorentzian metric with direction close to ei, p, or q, has a light 
cone emerging from B n _i(a) that includes X. This means that any time-like trajectory 
(geodesic) from apoint in X crosses B n _i(a) for all metrics with direction close to ei, 
p or q. See Fig. M where the light cone for p is shown to strictly include X. 
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Figure 6: Extended Geometry where complex geometric solutions are constructed. 

Now consider the slab to < x\ < t\. We prove a result on that slab and show that the 
Cauchy data at t\ are controlled so that the same estimate may be used on t% < X\ < t 2 
and on all of (0, 1) by induction. Let a± and fli be the two angles in (0, 2tt) such that 

cos ai% ) + sin a 1 9 ± (y ) = p, cos 0i9(y o ) + sin P^iVo) = q, 



where 9(x) 6 §™ is defined in (31) 



The complex geometric optics solutions are constructed as follows. We define v 
c =se p ' x and w = §le p ' x harmonic functions. Then we find that 



Vv = e k x \k\9(x), Vw = e k x \k\0 L (x) 



so that for the two harmonic functions v\ = cos a\V + sin a±w and w\ = cos piv+sm piw, 
we have on the slab < x\ < t\ that 



Vfi = cos a\9{x) + sin ai9 ± (x) = p + O^ 
V^i = cos f3 1 9(x) + sinf3 1 9 ± (x) = q + o(^ 



\k\ 

iV 
\k\ 

iV 



For t\ — i such that \k\t\ = ^ is sufficiently small, Vfi and Vwi for all x such that 
< X\ < t\ are two vector fields such that the associated Lorentzian metrics t)^p and 
b^~~ have ei as a time-like vector. 

Let us now assume that a is arbitrary but smooth. The main idea of CGO solutions 
is that we can construct solutions for arbitrary a that are close to the solutions corre- 
sponding to a = 1 for \k\ sufficiently large. We construct CGO solutions u p of Q (and 
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u p by replacing a by a) such that 



w p = ^e^(l + Vp), 

' (7 



with | & | "0,0 bounded in the C 1 norm since <r is sufficiently smooth by hypothesis. This 
result is proved in [9] following earlier work in [TT] . These solutions are constructed on 
M. n and then restricted to X; their boundary condition f p is therefore specified by the 
construction. For such a solution, we find that 

Vu p = ^=eP- x \p\(p + ct> p ), 



where |fc||0 p | is also bounded in the uniform norm. This shows that 

VS'Up(x) = 9(x) + <fi Pt i, V3fai p (x) = 9(x) + (f) P: r 

with | k 1 1 (p P} i | and |/c||0 Pir | bounded in the uniform norm. As a consequence, we have 
constructed solutions of ([I]) with a gradient that is close to the prescribed 9(x) corre- 
sponding to harmonic functions. Construct now the two linear combinations 

Vi tP = cosaifp + sinaiWp, Wi tP = cosf3iV p + sin/^u^, v p := Qu p , w p :=$t.Up. (34) 

Knowledge of the Cauchy data for v\ yP and w\ yP is inherited from that for v p and w p . 
Define vi )P and u>i iP similarly with a replaced by a. We choose \k\ sufficiently large and 
then t\\k\ sufficiently small so that (f) p is a negligible vector that does not perturb the 
Lorentzian metric much and so that 

Vv 1(P = p + 0(|fc|ti) + 0(M\k\~ 1 ) and Vw hp = q + 0(\k\tx) + 0{M\k\' 1 ) (35) 

are directions of Lorentzian metrics for which (i) ei is a time-like vector; and (ii) the 
light cone emerging from 5 n _i(a) includes X. Here, M is the uniform bound of a in 
Ff +3+E (r) E]. Note that this means that t\ should be chosen on the order of 
M\k\~ 2 once |fc| has been chosen so that M|/c| _1 is sufficiently small. 

The same properties hold for the vectors constructed by replacing a by a. Thus, the 



metric in (13) is given with a and /3 close to 1 and e(x) close to p for the function 
t>i iP and close to q for the function w ltP . Using Cauchy data on S := {x\ = 0}, we 
can then solve the linear equations on the slab 0\ := {0 = to < %\ < ti} and get the 
solution at the surface Si := {x\ = t\}. For the solutions V\ >p and Wx iP , we obtain as a 



slight modification of (20) [29J the stability result 



(/ \f-f\ 2 +\j-j\ 2 da+ [ E(dv)dx<c( [ \f-f\z + \j-j\*da+ [ \V5H\ 2 dx). 

v JEi J Ox V ~/So J Ox J 

(36) 

The above measurements 5H = H — H are those for the functions (v v i iP ) and 
(wi iP ,wi )P ). Such measurements can be constructed from the three measurements for 
v p , w p and v p + w p . This is the place where we use the three measurements stated in the 
theorem: we need to ensure that cr(x)|/iVf p + ^Vw p | 2 is available for any possible linear 
combination (/i, v) since the values of a\ and 0i will vary (and will be called ctj and 
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on each slab U < Xi < U+i. Since the measurements H for the 0— Laplacian problem 
are quadratic in the elliptic solution, three measurements are sufficient by polarization 
to allow us to construct a(x) |Vt>i >p | 2 and cr(x)\ Vu?i iP | 2 



On Ei, we have control on the Cauchy data of Vi :P and wi iP and hence of v p = Qu p 



and w p = tftiip thanks to (36) and (34). Here, we need that p and q be not too close 
to one- another (this is guaranteed by w < 1) so that the inversion of the 2x2 system 
is well-conditioned. On each slab, we define the angles ai and in order again to have 
Lorentzian metrics with directions close to p and q. We then obtain a similar estimate 



to (36) and pursue by induction until we reach the slab On '■= {tjv-i < Xi < t N = 1}. 

The stability results then apply to ^su p and $lu p , and we thus obtain a global estimate 
for a as in earlier sections. So far, the illuminations / prescribed on X to solve the elliptic 
problem are of a very specific type. In order for Qu p and 9?w p to be the solutions to 
the elliptic problems on X, then (fi, f 2 ) needs to be the trace of ($su p , $tu p ) on dX. 
It is for these illuminations that the three measurements Hij(x) for G / generate 
Lorentzian metrics that satisfy the above sufficient properties. Since a is not known, 
these traces are not known either. 

However, any Lorentzian metric that is sufficiently close to the Lorentzian metrics 
constructed with the real and imaginary parts of u p will inherit the same light cone prop- 
erties and, in particular, the fact that ei is a time-like vector for these new Lorentzian 
metrics throughout X = (0, 1) x B n _i(a). Therefore, there is an open set of boundary 



conditions (/i, f 2 ) close to (^su p \gx, $tu p \gx) such that the conclusion (36) holds, as well 



as the same expressions on the other slabs (9;. This concludes the proof of the result. □ 



Remark 4.6 The "three" measurements Hij for 6 I in (32) actually correspond 
to "two" physical measurements. Indeed, we can replace u e by U\- e and w_ e by u<i-- e in 
(J3]) and obtain in the limit crV«i • Vw 2 , which combined with H u and H 22 yields Hi 2 



defined in (32). The experimental acquisition of Hu is in fact sufficient to also acquire 



erVtti ■ Vu 2 as demonstrated in 

Remark 4.7 The above theorem shows a uniqueness and stability result for arbitrary, 
sufficiently smooth, conductivities. However, the boundary conditions / are quite spe- 
cific since they need to be sufficiently close to non-explicit, a-dependent, traces of com- 
plex geometrical optics solutions. In some sense, the difficulty inherent to the spatially 



varying Lorentzian metric f)(x) in (15) has been shifted to the difficulty of constructing 
adapted boundary conditions (illuminations). 

Note that the condition of flatness of the surfaces E« in the above construction is 
not essential. Surfaces with a geometry such as that depicted in Fig. [T] may also be 
considered. Such surfaces allow us to reduce the size of the domain X on which the 
conductivity o = 1 needs to be extended. Unless the domain X has a specific geometry 
similar to that of the domain O between Ei and E 2 in Fig. [TJ it seems necessary to 
augment the size of X to that of X as described above to obtain a global uniqueness 
result. 
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